home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1993 July / InfoMagic USENET CD-ROM July 1993.ISO / sources / misc / volume15 / ggems / part04 < prev    next >
Encoding:
Text File  |  1990-10-14  |  52.1 KB  |  1,898 lines

  1. Newsgroups: comp.sources.misc
  2. X-UNIX-From: craig@weedeater.math.yale.edu
  3. subject: v15i026: Graphics Gems, Part 4/5
  4. from: Craig Kolb <craig@weedeater.math.yale.edu>
  5. Sender: allbery@uunet.UU.NET (Brandon S. Allbery - comp.sources.misc)
  6.  
  7. Posting-number: Volume 15, Issue 26
  8. Submitted-by: Craig Kolb <craig@weedeater.math.yale.edu>
  9. Archive-name: ggems/part04
  10.  
  11. #! /bin/sh
  12. # This is a shell archive.  Remove anything before this line, then unpack
  13. # it by saving it into a file and typing "sh file".  To overwrite existing
  14. # files, type "sh file -c".  You can also feed this as standard input via
  15. # unshar, or by typing "sh <file", e.g..  If this archive is complete, you
  16. # will see the following message at the end:
  17. #        "End of archive 4 (of 5)."
  18. # Contents:  2DClip/cross.c AALines/AALines.c AALines/AATables.c
  19. #   AAPolyScan.c ConcaveScan.c Forms.c Interleave.c Sturm/sturm.c
  20. # Wrapped by craig@weedeater on Fri Oct 12 15:53:14 1990
  21. PATH=/bin:/usr/bin:/usr/ucb ; export PATH
  22. if test -f '2DClip/cross.c' -a "${1}" != "-c" ; then 
  23.   echo shar: Will not clobber existing file \"'2DClip/cross.c'\"
  24. else
  25. echo shar: Extracting \"'2DClip/cross.c'\" \(6788 characters\)
  26. sed "s/^X//" >'2DClip/cross.c' <<'END_OF_FILE'
  27. X
  28. X/*
  29. X * file cross.c:
  30. X *    calculate the intersections
  31. X */
  32. X#include    <math.h>
  33. X#include    "GraphicsGems.h"
  34. X#include    "line.h"
  35. X#include    "box.h"
  36. X/*
  37. X * cross_calc:
  38. X *
  39. X *    PURPOSE
  40. X *    calculate the intersections between the polygon
  41. X *    stored in poly and the linesegment stored in l
  42. X *    and put the intersections into psol.
  43. X *
  44. X *    poly    pointer to the structure containing the polygon
  45. X *    l    pointer to the structure containing the linesegment
  46. X *    psol    pointer to the pointer where intersections are stored
  47. X *    nsol    current number of intersections stored
  48. X *    nsmax    maximum storage in psol for intersections
  49. X *        if nsol exceeds nsmax additional storage is allocated
  50. X *
  51. X */
  52. Xcross_calc(poly, l, psol, nsol, nsmax)
  53. XCONTOUR    *poly;
  54. XSEGMENT    *l;
  55. XCLIST    **psol;
  56. Xshort    *nsol, nsmax;
  57. X{
  58. X    SEGMENT    *p;
  59. X    CLIST    *sol;
  60. X    double    s;
  61. X    long    x, y, a, b, c;
  62. X    int    psort(), type;
  63. X
  64. X    sol = *psol;
  65. X    p = poly->_s;
  66. X    do {
  67. X/*
  68. X * calculate the a, b and c coefficients and determine the
  69. X * type of intersection
  70. X */
  71. X
  72. X        a = (p->_to._y - p->_from._y)*(l->_to._x - l->_from._x) -
  73. X            (p->_to._x - p->_from._x)*(l->_to._y - l->_from._y);
  74. X        b = (p->_from._x - l->_from._x)*(l->_to._y - l->_from._y) -
  75. X            (p->_from._y - l->_from._y)*(l->_to._x - l->_from._x);
  76. X        c = (p->_from._x - l->_from._x)*(p->_to._y - p->_from._y) -
  77. X            (p->_from._y - l->_from._y)*(p->_to._x - p->_from._x);
  78. X        if(a == 0)
  79. X            type = (b == 0) ? COINCIDE : PARALLEL;
  80. X        else {
  81. X            if(a > 0) {
  82. X                if((b >= 0 && b <= a) &&
  83. X                   (c >= 0 && c <= a))
  84. X                    type = CROSS;
  85. X                else
  86. X                    type = NO_CROSS;
  87. X            }
  88. X            else {
  89. X                if((b <= 0 && b >= a) &&
  90. X                   (c <= 0 && c >= a))
  91. X                    type = CROSS;
  92. X                else
  93. X                    type = NO_CROSS;
  94. X            }
  95. X        }
  96. X/*
  97. X * process the interscetion found
  98. X */
  99. X        switch(type) {
  100. X            case NO_CROSS: case PARALLEL:
  101. X                break;
  102. X
  103. X            case CROSS:
  104. X                if(b == a || c == a || c == 0)
  105. X                    break;
  106. X                if(b == 0 && 
  107. X                   p_where(&(p->_prev->_from), &(p->_to), l) >= 0)
  108. X                    break;
  109. X                s = (double)b/(double)a;
  110. X                if(l->_from._x == l->_to._x)
  111. X                    x = l->_from._x;
  112. X                else
  113. X                    x = p->_from._x + 
  114. X                       (int)((p->_to._x - p->_from._x)*s);
  115. X                if(l->_from._y == l->_to._y)
  116. X                    y = l->_from._y;
  117. X                else
  118. X                    y = p->_from._y + 
  119. X                       (int)((p->_to._y - p->_from._y)*s);
  120. X
  121. X                if(*nsol == nsmax) {
  122. X                    nsmax *= 2;
  123. X                    *psol = sol = (CLIST *) realloc(sol,                             nsmax*sizeof(CLIST)); 
  124. X                }
  125. X                sol[*nsol]._p._x = x;
  126. X                sol[*nsol]._p._y = y;
  127. X                sol[*nsol]._type = STD;
  128. X                *nsol += 1;
  129. X                break;
  130. X
  131. X            case COINCIDE:
  132. X                if(p_where(&(p->_prev->_from),
  133. X                     &(p->_next->_to), l) > 0)
  134. X                    break;
  135. X                if(l->_from._x != l->_to._x) {
  136. X                    if((max(l->_from._x, l->_to._x) <
  137. X                        min(p->_from._x, p->_to._x)  ) ||
  138. X                       (min(l->_from._x, l->_to._x) >
  139. X                        max(p->_from._x, p->_to._x))   )
  140. X                        break;
  141. X                    if(min(l->_from._x, l->_to._x) <
  142. X                       min(p->_from._x, p->_to._x) ) {
  143. X                        if(*nsol == nsmax) {
  144. X                            nsmax *= 2;
  145. X                            *psol = sol = (CLIST *) realloc(sol,                                 nsmax*sizeof(CLIST));
  146. X                        }
  147. X                        sol[*nsol]._p._x = p->_from._x;
  148. X                        sol[*nsol]._p._y = p->_from._y;
  149. X                        sol[*nsol]._type = DELAY;
  150. X                        *nsol += 1;
  151. X                    }
  152. X                    if(max(l->_from._x, l->_to._x) >
  153. X                       max(p->_from._x, p->_to._x) ) {
  154. X                        if(*nsol == nsmax) {
  155. X                            nsmax *= 2;
  156. X                            *psol = sol = (CLIST *) realloc(sol,                                 nsmax*sizeof(CLIST));
  157. X                        }
  158. X                        sol[*nsol]._p._x = p->_to._x;
  159. X                        sol[*nsol]._p._y = p->_to._y;
  160. X                        sol[*nsol]._type = DELAY;
  161. X                        *nsol += 1;
  162. X                    }
  163. X                }
  164. X                else {
  165. X
  166. X                    if((max(l->_from._y, l->_to._y) <
  167. X                        min(p->_from._y, p->_to._y)  ) ||
  168. X                       (min(l->_from._y, l->_to._y) >
  169. X                        max(p->_from._y, p->_to._y)) )
  170. X                        break;
  171. X                    if(min(l->_from._y, l->_to._y) <
  172. X                       min(p->_from._y, p->_to._y) ) {
  173. X                        if(*nsol == nsmax) {
  174. X                            nsmax *= 2;
  175. X                            *psol = sol = (CLIST *) realloc(sol,                                 nsmax*sizeof(CLIST));
  176. X                        }
  177. X                        sol[*nsol]._p._x = p->_from._x;
  178. X                        sol[*nsol]._p._y = p->_from._y;
  179. X                        sol[*nsol]._type = DELAY;
  180. X                        *nsol += 1;
  181. X                    }
  182. X                    if(max(l->_from._y, l->_to._y) >
  183. X                       max(p->_from._y, p->_to._y) ) {
  184. X                        if(*nsol == nsmax) {
  185. X                            nsmax *= 2;
  186. X                            *psol = sol = (CLIST *) realloc(sol,                                 nsmax*sizeof(CLIST));
  187. X                        }
  188. X                        sol[*nsol]._p._x = p->_to._x;
  189. X                        sol[*nsol]._p._y = p->_to._y;
  190. X                        sol[*nsol]._type = DELAY;
  191. X                        *nsol += 1;
  192. X                    }
  193. X                }
  194. X                break;
  195. X        }
  196. X        p = p->_next;
  197. X    } while(p != poly->_s);
  198. X    qsort(sol, *nsol, sizeof(CLIST), psort);
  199. X}
  200. X
  201. X
  202. X/*
  203. X * p_where
  204. X *
  205. X *    PURPOSE
  206. X *    determine position of point p1 and p2 relative to
  207. X *    linesegment l. 
  208. X *    Return value
  209. X *    < 0    p1 and p2 lie at different sides of l
  210. X *    = 0    one of both points lie on l
  211. X *    > 0    p1 and p2 lie at same side of l
  212. X *
  213. X *    p1    pointer to coordinates of point
  214. X *    p2    pointer to coordinates of point
  215. X *    l    pointer to linesegment
  216. X * 
  217. X */
  218. Xp_where(p1, p2, l)
  219. XPOINT    *p1, *p2;
  220. XSEGMENT    *l;
  221. X{
  222. X    long    dx, dy, dx1, dx2, dy1, dy2, p_1, p_2;
  223. X
  224. X    dx  = l->_to._x - l->_from._x;
  225. X    dy  = l->_to._y - l->_from._y;
  226. X    dx1 = p1->_x - l->_from._x;
  227. X    dy1 = p1->_y - l->_from._y;
  228. X    dx2 = p2->_x - l->_to._x;
  229. X    dy2 = p2->_y - l->_to._y;
  230. X    p_1 = (dx*dy1 - dy*dx1);
  231. X    p_2 = (dx*dy2 - dy*dx2);
  232. X    if(p_1 == 0 || p_2 == 0)
  233. X        return(0);
  234. X    else {
  235. X        if((p_1 > 0 && p_2 < 0) || (p_1 < 0 && p_2 > 0))
  236. X            return(-1);
  237. X        else
  238. X            return(1);
  239. X    }
  240. X}
  241. X
  242. X
  243. X/*
  244. X * p_inside
  245. X *
  246. X *    PURPOSE
  247. X *    determine if the point stored in pt lies inside
  248. X *    the polygon stored in p
  249. X *    Return value:
  250. X *    FALSE    pt lies outside p
  251. X *    TRUE    pt lies inside  p
  252. X *
  253. X *    p    pointer to the polygon
  254. X *    pt    pointer to the point
  255. X */    
  256. Xboolean    p_inside(p, pt)
  257. XCONTOUR    *p;
  258. XPOINT    *pt;
  259. X{
  260. X    SEGMENT    l, *sp;
  261. X    CLIST    *sol;
  262. X    short    nsol = 0, nsmax = 2, reduce = 0, i;
  263. X    boolean    on_contour(), odd;
  264. X    
  265. X    l._from._x = p->_minx-2;
  266. X    l._from._y = pt->_y;
  267. X    l._to._x   = pt->_x;
  268. X    l._to._y   = pt->_y;
  269. X    sol = (CLIST *) calloc(2, sizeof(CLIST));
  270. X    cross_calc(p, &l, &sol, &nsol, nsmax);
  271. X    for(i=0; i<nsol-1; i++)
  272. X        if(sol[i]._type == DELAY && sol[i+1]._type == DELAY)
  273. X            reduce++;
  274. X    free(sol);
  275. X    odd = (nsol - reduce) & 0x01;
  276. X    return(odd ? !on_contour(p, pt) : FALSE);
  277. X}
  278. X
  279. X/*
  280. X * function used for sorting
  281. X */
  282. Xpsort(p1, p2)
  283. XCLIST    *p1, *p2;
  284. X{
  285. X    if(p1->_p._x != p2->_p._x)
  286. X        return(p1->_p._x - p2->_p._x);
  287. X    else 
  288. X        return(p1->_p._y - p2->_p._y);
  289. X}
  290. X
  291. X
  292. X/*
  293. X * on_contour
  294. X *
  295. X *    PURPOSE
  296. X *    determine if the point pt lies on the
  297. X *    contour p.
  298. X *    Return value
  299. X *    TRUE     point lies on contour
  300. X *    FALSE    point lies not on contour
  301. X *
  302. X *    p    pointer to the polygon structure
  303. X *    pt    pointer to the point
  304. X */
  305. Xboolean    on_contour(p, pt)
  306. XCONTOUR    *p;
  307. XPOINT    *pt;
  308. X{
  309. X    SEGMENT    *sp;
  310. X    long    dx1, dy1, dx2, dy2;
  311. X
  312. X    sp = p->_s;
  313. X    do {
  314. X        if((pt->_x >= min(sp->_from._x, sp->_to._x)) &&
  315. X           (pt->_x <= max(sp->_from._x, sp->_to._x))    ) { 
  316. X            dx1 = pt->_x - sp->_from._x;
  317. X            dx2 = sp->_to._x - pt->_x;
  318. X            dy1 = pt->_y - sp->_from._y;
  319. X            dy2 = sp->_to._y - pt->_y;
  320. X            if(dy1*dx2 == dy2*dx1)
  321. X                return(TRUE);
  322. X        }
  323. X        sp = sp->_next;
  324. X    } while(sp != p->_s);
  325. X    return(FALSE);
  326. X}
  327. END_OF_FILE
  328. if test 6788 -ne `wc -c <'2DClip/cross.c'`; then
  329.     echo shar: \"'2DClip/cross.c'\" unpacked with wrong size!
  330. fi
  331. # end of '2DClip/cross.c'
  332. fi
  333. if test -f 'AALines/AALines.c' -a "${1}" != "-c" ; then 
  334.   echo shar: Will not clobber existing file \"'AALines/AALines.c'\"
  335. else
  336. echo shar: Extracting \"'AALines/AALines.c'\" \(5139 characters\)
  337. sed "s/^X//" >'AALines/AALines.c' <<'END_OF_FILE'
  338. X/*  FILENAME: AALines.c  [revised 17 AUG 90]
  339. X
  340. X    AUTHOR:  Kelvin Thompson
  341. X
  342. X    DESCRIPTION:  Code to render an anti-aliased line, from
  343. X      "Rendering Anti-Aliased Lines" in _Graphics_Gems_.
  344. X
  345. X      This is almost exactly the code printed on pages 690-693
  346. X      of _Graphics_Gems_.  An overview of the code is on pages
  347. X      105-106.
  348. X
  349. X    LINK WITH:
  350. X      AALines.h -- Shared tables, symbols, etc.
  351. X      AAMain.c -- Calling code for this subroutine.
  352. X      AATables.c -- Initialize lookup tables.
  353. X*/
  354. X
  355. X
  356. X#include "AALines.h"
  357. X
  358. X#define SWAPVARS(v1,v2) ( temp=v1, v1=v2, v2=temp )
  359. X
  360. X#define FIXMUL(f1,f2) \
  361. X  ( \
  362. X  (((f1&0x0000ffff) * (f2&0x0000ffff)) >> 16) + \
  363. X  ((f1&0xffff0000)>>16) * (f2&0x0000ffff) + \
  364. X  ((f2&0xffff0000)>>16) * (f1&0x0000ffff) + \
  365. X  ((f1&0xffff0000)>>16) * (f2&0xffff0000) \
  366. X  )
  367. X
  368. X/* HARDWARE ASSUMPTIONS:
  369. X/*    * 32-bit, signed ints
  370. X/*    * 8-bit pixels, with initialized color table
  371. X/*    * pixels are memory mapped in a rectangular fashion */
  372. X
  373. X/* FIXED-POINT DATA TYPE */
  374. X#ifndef FX_FRACBITS
  375. X  typedef int FX;
  376. X# define FX_FRACBITS 16  /* bits of fraction in FX format */
  377. X# define FX_0        0   /* zero in fixed-point format */
  378. X#endif
  379. X
  380. X/* ASSUMED MACROS:
  381. X/*   SWAPVARS(v1,v2) -- swaps the contents of two variables
  382. X/*   PIXADDR(x,y) -- returns address of pixel at (x,y)
  383. X/*   COVERAGE(FXdist) -- lookup macro for pixel coverage 
  384. X        given perpendicular distance; takes a fixed-point
  385. X        integer and returns an integer in the range [0,255]
  386. X/*   SQRTFUNC(FXval) -- lookup macro for sqrt(1/(1+FXval^2))
  387. X        accepts and returns fixed-point numbers
  388. X/*   FIXMUL(FX1,FX2) -- multiplies two fixed-point numbers
  389. X        and returns the product as a fixed-point number   */
  390. X
  391. X/* BLENDING FUNCTION:
  392. X/*  'cover' is coverage -- in the range [0,255]
  393. X/*  'back' is background color -- in the range [0,255] */
  394. X#define BLEND(cover,back) ((((255-(cover))*(back))>>8)+(cover))
  395. X
  396. X/* LINE DIRECTION bits and tables */
  397. X#define DIR_STEEP  1  /* set when abs(dy) > abs(dx) */
  398. X#define DIR_NEGY   2  /* set whey dy < 0 */
  399. X
  400. X
  401. X/* pixel increment values 
  402. X/*   -- assume PIXINC(dx,dy) is a macro such that:
  403. X/*   PIXADDR(x0,y0) + PIXINC(dx,dy) = PIXADDR(x0+dx,y0+dy)  */
  404. Xstatic int adj_pixinc[4] = 
  405. X      { PIXINC(1,0), PIXINC(0,1), PIXINC(1,0), PIXINC(0,-1) };
  406. Xstatic int diag_pixinc[4] = 
  407. X      { PIXINC(1,1), PIXINC(1,1), PIXINC(1,-1), PIXINC(1,-1) };
  408. Xstatic int orth_pixinc[4] = 
  409. X      { PIXINC(0,1), PIXINC(1,0), PIXINC(0,-1), PIXINC(1,0) };
  410. X
  411. X/* Global 'Pmax' is initialized elsewhere.  It is the
  412. X   "maximum perpendicular distance" -- the sum of half the
  413. X   line width and the effective pixel radius -- in fixed format */
  414. XFX Pmax;
  415. X
  416. X
  417. X/***************  FUNCTION ANTI_LINE  ***************/
  418. X
  419. Xvoid Anti_Line ( X1, Y1, X2, Y2 )
  420. Xint X1, Y1, X2, Y2;
  421. X{
  422. Xint     Bvar,     /* decision variable for Bresenham's */
  423. X        Bainc,   /* adjacent-increment for 'Bvar' */
  424. X        Bdinc;   /* diagonal-increment for 'Bvar' */
  425. XFX         Pmid,      /* perp distance at Bresenham's pixel */
  426. X           Pnow,      /* perp distance at current pixel (ortho loop) */
  427. X           Painc,     /* adjacent-increment for 'Pmid' */
  428. X           Pdinc,     /* diagonal-increment for 'Pmid' */
  429. X           Poinc;     /* orthogonal-increment for 'Pnow'--also equals 'k' */
  430. Xchar     *mid_addr,   /* pixel address for Bresenham's pixel */
  431. X         *now_addr;   /* pixel address for current pixel */
  432. Xint     addr_ainc,   /* adjacent pixel address offset */
  433. X        addr_dinc,   /* diagonal pixel address offset */
  434. X        addr_oinc;   /* orthogonal pixel address offset */
  435. Xint dx,dy,dir;        /* direction and deltas */
  436. XFX slope;            /* slope of line */
  437. Xint temp;
  438. X
  439. X/* rearrange ordering to force left-to-right */
  440. Xif     ( X1 > X2 )
  441. X      { SWAPVARS(X1,X2);  SWAPVARS(Y1,Y2); }
  442. X
  443. X/* init deltas */
  444. Xdx = X2 - X1;  /* guaranteed non-negative */
  445. Xdy = Y2 - Y1;
  446. X
  447. X
  448. X/* calculate direction (slope category) */
  449. Xdir = 0;
  450. Xif ( dy < 0 )   { dir |= DIR_NEGY;  dy = -dy; }
  451. Xif ( dy > dx )  { dir |= DIR_STEEP; SWAPVARS(dx,dy); }
  452. X
  453. X/* init address stuff */
  454. Xmid_addr = PIXADDR(X1,Y1);
  455. Xaddr_ainc = adj_pixinc[dir];
  456. Xaddr_dinc = diag_pixinc[dir];
  457. Xaddr_oinc = orth_pixinc[dir];
  458. X
  459. X/* perpendicular measures */
  460. Xslope =  (dy << FX_FRACBITS) / dx;
  461. XPoinc = SQRTFUNC( slope );
  462. XPainc = FIXMUL( slope, Poinc );
  463. XPdinc = Painc - Poinc;
  464. XPmid = FX_0;
  465. X
  466. X/* init Bresenham's */
  467. XBainc = dy << 1;
  468. XBdinc = (dy-dx) << 1;
  469. XBvar = Bainc - dx;
  470. X
  471. Xdo
  472. X      {
  473. X      /* do middle pixel */
  474. X      *mid_addr = BLEND( COVERAGE(abs(Pmid)), *mid_addr );
  475. X
  476. X      /* go up orthogonally */
  477. X      for (
  478. X          Pnow = Poinc-Pmid,  now_addr = mid_addr+addr_oinc;
  479. X          Pnow < Pmax;
  480. X          Pnow += Poinc,      now_addr += addr_oinc
  481. X          )
  482. X       *now_addr = BLEND( COVERAGE(Pnow), *now_addr );
  483. X
  484. X      /* go down orthogonally */
  485. X      for (
  486. X          Pnow = Poinc+Pmid,  now_addr = mid_addr-addr_oinc;
  487. X          Pnow < Pmax;
  488. X          Pnow += Poinc,      now_addr -= addr_oinc
  489. X          )
  490. X       *now_addr = BLEND( COVERAGE(Pnow), *now_addr );
  491. X
  492. X
  493. X      /* update Bresenham's */
  494. X      if ( Bvar < 0 )
  495. X        {
  496. X        Bvar += Bainc;
  497. X        mid_addr = (char *) ((int)mid_addr + addr_ainc);
  498. X        Pmid += Painc;
  499. X        }
  500. X      else
  501. X        {
  502. X        Bvar += Bdinc;
  503. X        mid_addr = (char *) ((int)mid_addr + addr_dinc);
  504. X        Pmid += Pdinc;
  505. X        }
  506. X
  507. X      --dx;
  508. X      } while ( dx >= 0 );
  509. X}
  510. END_OF_FILE
  511. if test 5139 -ne `wc -c <'AALines/AALines.c'`; then
  512.     echo shar: \"'AALines/AALines.c'\" unpacked with wrong size!
  513. fi
  514. # end of 'AALines/AALines.c'
  515. fi
  516. if test -f 'AALines/AATables.c' -a "${1}" != "-c" ; then 
  517.   echo shar: Will not clobber existing file \"'AALines/AATables.c'\"
  518. else
  519. echo shar: Extracting \"'AALines/AATables.c'\" \(5772 characters\)
  520. sed "s/^X//" >'AALines/AATables.c' <<'END_OF_FILE'
  521. X/*  FILENAME: AATables.c  [revised 18 AUG 90]
  522. X
  523. X    DESCRIPTION:  Initialization of lookup tables and frame buffer
  524. X      for anti-aliased line rendering demo.
  525. X
  526. X    LINK WITH:
  527. X      AALines.h -- Shared variables and symbols for renderer.
  528. X      AAMain.c -- Calling routine for renderer.
  529. X      AALines.c -- Anti-aliased line rendering code.
  530. X*/
  531. X
  532. X#include <math.h>
  533. X#include "AALines.h"
  534. X
  535. X/* programs in this file */
  536. Xextern void Anti_Init();
  537. Xstatic void Sqrt_Init();
  538. X
  539. X/* globals defined here */
  540. Xchar *fbuff;
  541. XUFX *sqrtfunc=0;
  542. Xint sqrtcells=1024;
  543. Xint sqrtshift;
  544. X
  545. X/* AA sizes */
  546. Xfloat line_r=1.0;     /* line radius */
  547. Xfloat pix_r=SQRT_2;   /* pixel radius */
  548. XFX *coverage=0;
  549. Xint covercells=128;
  550. Xint covershift;
  551. X
  552. X
  553. X
  554. X
  555. X/*  ************************************
  556. X    *                                  *
  557. X    *          Anti_Init               *
  558. X    *                                  *
  559. X    ************************************
  560. X
  561. X    DESCRIPTION:  Initialize everything for the anti-aliased line
  562. X      renderer in 'AALines.c':  allocate the frame buffer,
  563. X      set up lookup tables, etc.
  564. X
  565. X      For hints about initializing the coverage table, see 
  566. X      "Area of Intersection: Circle and a Thick Line" on pages
  567. X      40-42 of _Graphics_Gems_.
  568. X*/
  569. X
  570. X#define FLOAT_TO_CELL(flt)  ((int) ((flt) * 255.0 + 0.5))
  571. X#define MAXVAL_CELL         255
  572. X
  573. Xvoid Anti_Init ( )
  574. X{
  575. Xint *thiscell;
  576. Xdouble maxdist,nowdist,incdist;
  577. Xint tablebits,radbits;
  578. Xint tablecells;
  579. Xstatic int tablesize=0;
  580. Xdouble fnear,ffar,fcover;
  581. Xdouble half,invR,invpiRsq,invpi,Rsq;
  582. Xdouble sum_r;
  583. Xdouble inv_log_2;
  584. Xextern FX Pmax;
  585. X
  586. X/* alloc & init frame buffer */
  587. Xfbuff = (char *) malloc( xpix*ypix );
  588. X  {
  589. X  register int i;
  590. X  for ( i=xpix*ypix-1; i>=0; --i )  fbuff[i] = 0;
  591. X  }
  592. X
  593. X/* init */
  594. Xinv_log_2 = 1.0 / log( 2.0 );
  595. Xsum_r = line_r + pix_r;
  596. Xtablebits = (int) ( log((double)covercells) * inv_log_2 + 0.99 );
  597. Xradbits = (int) ( log((double)sum_r) * inv_log_2 ) + 1;
  598. Xcovershift = FX_FRACBITS - (tablebits-radbits);
  599. X
  600. X/* constants */
  601. Xhalf = 0.5;
  602. XinvR = 1.0 / pix_r;
  603. Xinvpi = 1.0 / PI;
  604. XinvpiRsq = invpi * invR * invR;
  605. XRsq = pix_r * pix_r;
  606. X#define FRACCOVER(d) (half - d*sqrt(Rsq-d*d)*invpiRsq - invpi*asin(d*invR))
  607. X
  608. X/* allocate table */
  609. XPmax = FLOAT_TO_FX(sum_r);
  610. XPmax >>= covershift;
  611. Xtablecells = Pmax + 2;
  612. XPmax <<= covershift;
  613. Xif ( coverage  &&  tablecells > tablesize )
  614. X  { free( coverage ); coverage = 0;  tablesize = 0; }
  615. Xif ( coverage == 0 )
  616. X  {
  617. X  coverage = (FX *) malloc( tablecells * sizeof(int) );
  618. X  tablesize = tablecells;
  619. X  }
  620. X
  621. X/* init for fill loops */
  622. Xnowdist = 0.0;
  623. Xthiscell = coverage;
  624. Xincdist = sum_r / (double)(tablecells-2);
  625. X
  626. X/* fill fat portion */
  627. Xif ( pix_r <= line_r )
  628. X  {
  629. X  maxdist = line_r - pix_r;
  630. X  for ( ;
  631. X      nowdist <= maxdist;
  632. X      nowdist += incdist,  ++thiscell
  633. X      )
  634. X    {
  635. X    *thiscell = MAXVAL_CELL;
  636. X    }
  637. X  }
  638. X
  639. X/* fill skinny portion */
  640. Xelse
  641. X  {
  642. X  /* loop till edge of line, or end of skinny, whichever comes first */
  643. X  maxdist = pix_r - line_r;
  644. X  if ( maxdist > line_r )
  645. X    maxdist = line_r;
  646. X  for ( ;
  647. X      nowdist < maxdist;
  648. X      nowdist += incdist,  ++thiscell
  649. X      )
  650. X    {
  651. X    fnear = line_r - nowdist;
  652. X    ffar = line_r + nowdist;
  653. X    fcover = 1.0 - FRACCOVER(fnear) - FRACCOVER(ffar);
  654. X    *thiscell = FLOAT_TO_CELL(fcover);
  655. X    }
  656. X
  657. X  /* loop till end of skinny -- only run on super-skinny */
  658. X  maxdist = pix_r - line_r;
  659. X  for ( ;
  660. X      nowdist < maxdist;
  661. X      nowdist += incdist,  ++thiscell
  662. X      )
  663. X    {
  664. X    fnear = nowdist - line_r;
  665. X    ffar = nowdist + line_r;
  666. X    fcover = FRACCOVER(fnear) - FRACCOVER(ffar);
  667. X    *thiscell = FLOAT_TO_CELL(fcover);
  668. X    }
  669. X  }
  670. X
  671. X/* loop till edge of line */
  672. Xmaxdist = line_r;
  673. Xfor ( ;
  674. X    nowdist < maxdist;
  675. X    nowdist += incdist,  ++thiscell
  676. X    )
  677. X  {
  678. X  fnear = line_r - nowdist;
  679. X  fcover = 1.0 - FRACCOVER(fnear);
  680. X  *thiscell = FLOAT_TO_CELL(fcover);
  681. X  }
  682. X
  683. X/* loop till max separation */
  684. Xmaxdist = line_r + pix_r;
  685. Xfor ( ;
  686. X    nowdist < maxdist;
  687. X    nowdist += incdist,  ++thiscell
  688. X    )
  689. X  {
  690. X  fnear = nowdist - line_r;
  691. X  fcover = FRACCOVER(fnear);
  692. X  *thiscell = FLOAT_TO_CELL(fcover);
  693. X  }
  694. X
  695. X/* finish off table */
  696. X*thiscell = FLOAT_TO_CELL(0.0);
  697. Xcoverage[tablecells-1] = FLOAT_TO_CELL(0.0);
  698. X
  699. XSqrt_Init();
  700. X}
  701. X
  702. X
  703. X
  704. X
  705. X/*  *******************************
  706. X    *                             *
  707. X    *       Sqrt_Init             *
  708. X    *                             *
  709. X    *******************************
  710. X
  711. X    DESCRIPTION:  Initialize the lookup table for the function
  712. X      sqrt(1/(1+x^2))).  The table takes a shifted fixed-point
  713. X      value as an index and returns a fixed-point value.  Input
  714. X      values are in the range [0,1] inclusive.  The number of
  715. X      cells in the table is a power of two plus one (the extra
  716. X      cell provides an entry for an input value of exactly 1).
  717. X
  718. X    GLOBALS:
  719. X      sqrtcells -- Number of cells to use in the table (must
  720. X    be set before calling this routine).  This number is
  721. X    rounded up to the nearest power of two (the global
  722. X    variable itself is unchanged).
  723. X      sqrtshift -- Bits to shift a fixed-point (FX) number
  724. X    to generate a table index.
  725. X      sqrtfunc -- Lookup table for the function.
  726. X*/
  727. X
  728. Xstatic void Sqrt_Init ( )
  729. X{
  730. XUFX *thiscell;
  731. Xdouble nowval,incval;
  732. Xint tablebits;
  733. Xint tablecells;
  734. Xdouble one;
  735. X
  736. X/* init */
  737. Xtablebits = (int) ( log((double)sqrtcells) / log(2.0) + 0.999 );
  738. Xtablecells = (1 << tablebits) + 1;  /* one more that requested */
  739. Xsqrtshift = FX_FRACBITS - tablebits;
  740. Xone = 1.0;
  741. X
  742. X/* allocate table */
  743. Xif ( sqrtfunc )
  744. X  free( sqrtfunc );
  745. Xsqrtfunc = (UFX *) malloc( tablecells * sizeof(int) );
  746. X
  747. X/* init for fill loop */
  748. Xincval = one / (double)(tablecells-1);  /* a negative power of two */
  749. X
  750. Xfor (
  751. X    nowval = 0.0,      thiscell = sqrtfunc;
  752. X    nowval < 1.0;
  753. X    nowval += incval,  ++thiscell
  754. X    )
  755. X  {
  756. X  *thiscell = FLOAT_TO_FX( sqrt(one/(one+nowval*nowval)) );
  757. X  }
  758. X
  759. Xsqrtfunc[tablecells-1] = FLOAT_TO_FX( sqrt(0.5) );
  760. X}
  761. END_OF_FILE
  762. if test 5772 -ne `wc -c <'AALines/AATables.c'`; then
  763.     echo shar: \"'AALines/AATables.c'\" unpacked with wrong size!
  764. fi
  765. # end of 'AALines/AATables.c'
  766. fi
  767. if test -f 'AAPolyScan.c' -a "${1}" != "-c" ; then 
  768.   echo shar: Will not clobber existing file \"'AAPolyScan.c'\"
  769. else
  770. echo shar: Extracting \"'AAPolyScan.c'\" \(7519 characters\)
  771. sed "s/^X//" >'AAPolyScan.c' <<'END_OF_FILE'
  772. X/* 
  773. XFast Anti-Aliasing Polygon Scan Converstion
  774. Xby Jack Morrison
  775. Xfrom "Graphics Gems", Academic Press, 1990
  776. X*/
  777. X
  778. X/*
  779. X* Anti-aliased polygon scan conversion by Jack Morrison
  780. X*
  781. X* This code renders a polygon, computing subpixel coverage at
  782. X* 8 times Y and 16 times X display resolution for anti-aliasing.
  783. X* One optimization left out for clarity is the use of incremental
  784. X* interpolations. X coordinate interpolation in particular can be
  785. X* with integers. See Dan Field's article in ACM Transactions on
  786. X* Graphics, January 1985 for a fast incremental interpolator.
  787. X*/
  788. X#include "GraphicsGems.h"
  789. X
  790. X#define    SUBYRES    8        /* subpixel Y resolution per scanline */
  791. X#define    SUBXRES    16        /* subpixel X resolution per pixel */
  792. X#define    MAX_AREA    (SUBYRES*SUBXRES)
  793. X#define    MODRES(y)    ((y) & 7)        /*subpixel Y modulo */
  794. X#define MAX_X    0x7FFF    /* subpixel X beyond right edge */
  795. X
  796. Xtypedef struct SurfaceStruct {  /* object shading surface info */
  797. X    int    red, green, blue;           /* color components */
  798. X    } Surface;
  799. X/*
  800. X* In  real life, SurfaceStruct will contain many more parameters as
  801. X* required by the shading and rendering programs, such as diffuse
  802. X* and specular factors, texture information, transparency, etc.
  803. X*/
  804. X
  805. Xtypedef struct VertexStruct    {    /* polygon vertex */
  806. X    Vector3    model, world,        /* geometric information */
  807. X            normal, image;
  808. X    int y;                    /* subpixel display coordinate */
  809. X    } Vertex;
  810. X
  811. XVertex *Vleft, *VnextLeft;        /* current left edge */
  812. XVertex *Vright, *VnextRight;    /* current right edge */
  813. X
  814. Xstruct    SubPixel  {            /* subpixel extents for scanline */
  815. X    int xLeft, xRight;
  816. X    } sp[SUBYRES];
  817. X
  818. Xint    xLmin, xLmax;        /* subpixel x extremes for scanline */
  819. Xint    xRmax, xRmin;        /* (for optimization shortcut) */
  820. X
  821. X/* Compute sub-pixel x coordinate for vertex */
  822. Xextern int screenX(/* Vertex *v */);
  823. X
  824. X/* Interpolate vertex information */
  825. Xextern void vLerp(/* double alpha, Vertex *Va, *Vb, *Vout */);
  826. X
  827. X/* Render polygon for one pixel, given coverage area */
  828. X/*  and bitmask */
  829. Xextern void renderPixel(/* int x, y, Vertex *V,
  830. X                        int area, unsigned mask[], 
  831. X                        Surface *object */);
  832. X
  833. X/*
  834. X * Render shaded polygon
  835. X */
  836. XdrawPolygon(polygon, numVertex, object)
  837. X    Vertex    polygon[];        /*clockwise clipped vertex list */
  838. X    int    numVertex;            /*number of vertices in polygon */
  839. X
  840. X    Surface *object;            /* shading parms for this object */
  841. X{
  842. X    Vertex *endPoly;            /* end of polygon vertex list */
  843. X    Vertex VscanLeft, VscanRight;    /* interpolated vertices */                                 /* at scanline */
  844. X    double aLeft, aRight;            /* interpolation ratios */
  845. X    struct SubPixel *sp_ptr;        /* current subpixel info */
  846. X    int xLeft, xNextLeft;            /* subpixel coordinates for */
  847. X    int  xRight, xNextRight;        /* active polygon edges */
  848. X    int i,y;                        
  849. X
  850. X/* find vertex with minimum y (display coordinate) */
  851. XVleft = polygon;
  852. Xfor  (i=1; i<numVertex; i++)
  853. X    if  (polygon[i].y < Vleft->y)
  854. X        Vleft = &polygon[i];
  855. XendPoly = &polygon[numVertex-1];
  856. X
  857. X/* initialize scanning edges */
  858. XVright = VnextRight = VnextLeft = Vleft;
  859. X
  860. X/* prepare bottom of initial scanline - no coverage by polygon */
  861. Xfor (i=0; i<SUBYRES; i++)
  862. X    sp[i].xLeft = sp[i].xRight = -1;
  863. XxLmin = xRmin = MAX_X;
  864. XxLmax = xRmax = -1;
  865. X
  866. X/* scan convert for each subpixel from bottom to top */
  867. Xfor (y+Vleft->y; ; y++)    {
  868. X
  869. X    while (y == VnextLeft->y)    {    /* reached next left vertex */
  870. X        VnextLeft = (Vleft=VnextLeft) + 1;     /* advance */
  871. X        if (VnextLeft > endPoly)            /* (wraparound) */
  872. X            VnextLeft = polygon;
  873. X        if (VnextLeft == Vright)    /* all y's same?  */
  874. X            return;                /* (null polygon) */ 
  875. X        xLeft = screenX(Vleft);
  876. X        xNextLeft = screenX(VnextLeft);
  877. X    }
  878. X
  879. X    while (y == VnextRight->y)  { /*reached next right vertex */
  880. X        VnextRight = (Vright=VnextRight) -1;
  881. X        if (VnextRight < polygon)            /* (wraparound) */
  882. X            VnextRight = endPoly;
  883. X        xRight = screenX(Vright);
  884. X        xNextRight = screenX(VnextRight);
  885. X    }
  886. X
  887. X    if (y>VnextLeft->y || y>VnextRight->y)    {
  888. X                /* done, mark uncovered part of last scanline */
  889. X        for (; MODRES(y); y++)
  890. X            sp[MODRES(y)].xLeft = sp[MODRES(y)].xRight = -1;
  891. X        renderScanline(Vleft, Vright, y/SUBYRES, object);
  892. X        return;
  893. X    }
  894. X
  895. X/*
  896. X * Interpolate sub-pixel x endpoints at this y,
  897. X * and update extremes for pixel coherence optimization
  898. X */
  899. X    
  900. X    sp_ptr = &sp[MODRES(y)];
  901. X    aLeft = (double)(y - Vleft->y) / (VnextLeft->y - Vleft->y);
  902. X    sp_ptr->xLeft = LERP(aLeft, xLeft, xNextLeft);
  903. X    if (sp_ptr->xLeft < xLmin)
  904. X        xLmin = sp_ptr->xLeft;
  905. X    if (sp_ptr->xLeft > xLmax)
  906. X        xLmax = sp_ptr->xLeft;
  907. X
  908. X    aRight = (double)(y - Vright->y) / (VnextRight->y 
  909. X                    - Vright->y);
  910. X    sp_ptr->xRight = LERP(aRight, xRight, xNextRight);
  911. X    if (sp_ptr->xRight < xRmin)
  912. X        xRmin = sp_ptr->xRight;
  913. X    if (sp_ptr->xRight > xRmax)
  914. X        xRmax = sp_ptr->xRight;
  915. X
  916. X    if (MODRES(y) == SUBYRES-1)    {    /* end of scanline */
  917. X            /* interpolate edges to this scanline */
  918. X        vLerp(aLeft, Vleft, VnextLeft, &VscanLeft);
  919. X        vLerp(aRight, Vright, VnextRight, &VscanRight);
  920. X        renderScanline(&VscanLeft, &VscanRight, y/SUBYRES, object);
  921. X        xLmin = xRmin = MAX_X;         /* reset extremes */
  922. X        xLmax = xRmax = -1;
  923. X    }
  924. X  }
  925. X}
  926. X
  927. X/*
  928. X * Render one scanline of polygon
  929. X */
  930. X
  931. XrenderScanline(Vl, Vr, y, object)
  932. X    Vertex *Vl, *Vr;     /* polygon vertices interpolated */
  933. X                    /* at scanline */   
  934. X    int y;            /* scanline coordinate */
  935. X    Surface *object;    /* shading parms for this object */
  936. X{
  937. X    Vertex Vpixel;    /*object info interpolated at one pixel */
  938. X    unsigned mask[SUBYRES];    /*pixel coverage bitmask */
  939. X    int x;            /* leftmost subpixel of current pixel */
  940. X
  941. X    for (x=SUBXRES*FLOOR(xLmin/SUBXRES); x<=xRmax; x+=SUBXRES) {
  942. X        vLerp((double)(x-xLmin)/(xRmax-xLmin), Vl, Vr, &Vpixel);
  943. X        computePixelMask(x, mask);
  944. X        renderPixel(x/SUBXRES, y, &Vpixel,
  945. X                    /*computePixel*/Coverage(x), mask, object);
  946. X    }
  947. X}
  948. X
  949. X/*
  950. X * Compute number of subpixels covered by polygon at current pixel
  951. X*/
  952. X/*computePixel*/Coverage(x)
  953. X    int x;            /* left subpixel of pixel */
  954. X{
  955. X    int  area;            /* total covered area */
  956. X    int partialArea;      /* covered area for current subpixel y */
  957. X    int xr = x+SUBXRES-1;    /*right subpixel of pixel */
  958. X    int y;
  959. X
  960. X    /* shortcut for common case of fully covered pixel */
  961. X    if (x>xLmax && x<xRmin)
  962. X        return MAX_AREA;
  963. X    
  964. X    for (area=y=0; y<SUBYRES; y++) {
  965. X        partialArea = MIN(sp[y].xRight, xr)
  966. X             - MAX(sp[y].xLeft, x) + 1;
  967. X        if (partialArea > 0)
  968. X            area += partialArea;
  969. X    }
  970. X    return area;
  971. X}
  972. X
  973. X/* Compute bitmask indicating which subpixels are covered by 
  974. X * polygon at current pixel. (Not all hidden-surface methods
  975. X * need this mask. )
  976. X*/
  977. XcomputePixelMask(x, mask)
  978. X    int x;            /* left subpixel of pixel */
  979. X    unsigned mask[];    /* output bitmask */
  980. X{
  981. X    static unsigned leftMaskTable[] =
  982. X        { 0xFFFF, 0x7FFF, 0x3FFF, 0x1FFF, 0x0FFF, 0x07FF, 0x03FF,
  983. X          0x01FF, 0x00FF, 0x007F, 0x003F, 0x001F, 0x000F, 0x0007,
  984. X          0x0003, 0x0001  };
  985. X    static unsigned rightMaskTable[] = 
  986. X        { 0x8000, 0xC000, 0xE000, 0xF000, 0xF800, 0xFC00, 
  987. X          0xFE00, 0xFF00, 0xFF80, 0xFFC0, 0xFFE0, 0xFFF0,
  988. X          0xFFF8, 0xFFFC, 0xFFFE, 0xFFFF   };
  989. X    unsigned leftMask, rightMask;         /* partial masks */
  990. X    int xr = x+SUBXRES-1;            /* right subpixel of pixel */
  991. X    int y;
  992. X
  993. X/* shortcut for common case of fully covered pixel */
  994. X    if (x>xLmax && x<xRmin)     {
  995. X        for (y=0; y<SUBYRES; y++)
  996. X            mask[y] = 0xFFFF;
  997. X    } else     {
  998. X        for (y=0; y<SUBYRES; y++)    {
  999. X            if (sp[y].xLeft < x)     /* completely left of pixel*/
  1000. X                leftMask = 0xFFFF;
  1001. X            else if (sp[y].xLeft > xr)  /* completely right */    
  1002. X                leftMask = 0;
  1003. X            else
  1004. X                leftMask = leftMaskTable[sp[y].xLeft -x];
  1005. X
  1006. X            if (sp[y].xRight > xr)      /* completely  */
  1007. X                            /* right of pixel*/
  1008. X                rightMask = 0xFFFF;
  1009. X            else if (sp[y].xRight < x)    /*completely left */
  1010. X                rightMask = 0;
  1011. X            else
  1012. X                rightMask = rightMaskTable[sp[y].xRight -x];
  1013. X            mask[y] = leftMask & rightMask;
  1014. X        }
  1015. X    }
  1016. X}
  1017. END_OF_FILE
  1018. if test 7519 -ne `wc -c <'AAPolyScan.c'`; then
  1019.     echo shar: \"'AAPolyScan.c'\" unpacked with wrong size!
  1020. fi
  1021. # end of 'AAPolyScan.c'
  1022. fi
  1023. if test -f 'ConcaveScan.c' -a "${1}" != "-c" ; then 
  1024.   echo shar: Will not clobber existing file \"'ConcaveScan.c'\"
  1025. else
  1026. echo shar: Extracting \"'ConcaveScan.c'\" \(4931 characters\)
  1027. sed "s/^X//" >'ConcaveScan.c' <<'END_OF_FILE'
  1028. X/*
  1029. X * Concave Polygon Scan Conversion
  1030. X * by Paul Heckbert
  1031. X * from "Graphics Gems", Academic Press, 1990
  1032. X */
  1033. X
  1034. X/*
  1035. X * concave: scan convert nvert-sided concave non-simple polygon with vertices at
  1036. X * (point[i].x, point[i].y) for i in [0..nvert-1] within the window win by
  1037. X * calling spanproc for each visible span of pixels.
  1038. X * Polygon can be clockwise or counterclockwise.
  1039. X * Algorithm does uniform point sampling at pixel centers.
  1040. X * Inside-outside test done by Jordan's rule: a point is considered inside if
  1041. X * an emanating ray intersects the polygon an odd number of times.
  1042. X * drawproc should fill in pixels from xl to xr inclusive on scanline y,
  1043. X * e.g:
  1044. X *    drawproc(y, xl, xr)
  1045. X *    int y, xl, xr;
  1046. X *    {
  1047. X *        int x;
  1048. X *        for (x=xl; x<=xr; x++)
  1049. X *        pixel_write(x, y, pixelvalue);
  1050. X *    }
  1051. X *
  1052. X *  Paul Heckbert    30 June 81, 18 Dec 89
  1053. X */
  1054. X
  1055. X#include <stdio.h>
  1056. X#include <math.h>
  1057. X#include "GraphicsGems.h"
  1058. X
  1059. X#define ALLOC(ptr, type, n)  ASSERT(ptr = (type *)malloc((n)*sizeof(type)))
  1060. X
  1061. Xtypedef struct {        /* window: a discrete 2-D rectangle */
  1062. X    int x0, y0;            /* xmin and ymin */
  1063. X    int x1, y1;            /* xmax and ymax (inclusive) */
  1064. X} Window;
  1065. X
  1066. Xtypedef struct {        /* a polygon edge */
  1067. X    double x;    /* x coordinate of edge's intersection with current scanline */
  1068. X    double dx;    /* change in x with respect to y */
  1069. X    int i;    /* edge number: edge i goes from pt[i] to pt[i+1] */
  1070. X} Edge;
  1071. X
  1072. Xstatic int n;            /* number of vertices */
  1073. Xstatic Point2 *pt;        /* vertices */
  1074. X
  1075. Xstatic int nact;        /* number of active edges */
  1076. Xstatic Edge *active;        /* active edge list:edges crossing scanline y */
  1077. X
  1078. Xint compare_ind(), compare_active();
  1079. X
  1080. Xconcave(nvert, point, win, spanproc)
  1081. Xint nvert;            /* number of vertices */
  1082. XPoint2 *point;            /* vertices of polygon */
  1083. XWindow *win;            /* screen clipping window */
  1084. Xvoid (*spanproc)();        /* called for each span of pixels */
  1085. X{
  1086. X    int k, y0, y1, y, i, j, xl, xr;
  1087. X    int *ind;        /* list of vertex indices, sorted by pt[ind[j]].y */
  1088. X
  1089. X    n = nvert;
  1090. X    pt = point;
  1091. X    if (n<=0) return;
  1092. X    ALLOC(ind, int, n);
  1093. X    ALLOC(active, Edge, n);
  1094. X
  1095. X    /* create y-sorted array of indices ind[k] into vertex list */
  1096. X    for (k=0; k<n; k++)
  1097. X    ind[k] = k;
  1098. X    qsort(ind, n, sizeof ind[0], compare_ind);    /* sort ind by pt[ind[k]].y */
  1099. X
  1100. X    nact = 0;                /* start with empty active list */
  1101. X    k = 0;                /* ind[k] is next vertex to process */
  1102. X    y0 = MAX(win->y0, ceil(pt[ind[0]].y-.5));        /* ymin of polygon */
  1103. X    y1 = MIN(win->y1, floor(pt[ind[n-1]].y-.5));    /* ymax of polygon */
  1104. X
  1105. X    for (y=y0; y<=y1; y++) {        /* step through scanlines */
  1106. X    /* scanline y is at y+.5 in continuous coordinates */
  1107. X
  1108. X    /* check vertices between previous scanline and current one, if any */
  1109. X    for (; k<n && pt[ind[k]].y<=y+.5; k++) {
  1110. X        /* to simplify, if pt.y=y+.5, pretend it's above */
  1111. X        /* invariant: y-.5 < pt[i].y <= y+.5 */
  1112. X        i = ind[k];    
  1113. X        /*
  1114. X         * insert or delete edges before and after vertex i (i-1 to i,
  1115. X         * and i to i+1) from active list if they cross scanline y
  1116. X         */
  1117. X        j = i>0 ? i-1 : n-1;    /* vertex previous to i */
  1118. X        if (pt[j].y <= y-.5)    /* old edge, remove from active list */
  1119. X        delete(j);
  1120. X        else if (pt[j].y > y+.5)    /* new edge, add to active list */
  1121. X        insert(j, y);
  1122. X        j = i<n-1 ? i+1 : 0;    /* vertex next after i */
  1123. X        if (pt[j].y <= y-.5)    /* old edge, remove from active list */
  1124. X        delete(i);
  1125. X        else if (pt[j].y > y+.5)    /* new edge, add to active list */
  1126. X        insert(i, y);
  1127. X    }
  1128. X
  1129. X    /* sort active edge list by active[j].x */
  1130. X    qsort(active, nact, sizeof active[0], compare_active);
  1131. X
  1132. X    /* draw horizontal segments for scanline y */
  1133. X    for (j=0; j<nact; j+=2) {    /* draw horizontal segments */
  1134. X        /* span 'tween j & j+1 is inside, span tween j+1 & j+2 is outside */
  1135. X        xl = ceil(active[j].x-.5);        /* left end of span */
  1136. X        if (xl<win->x0) xl = win->x0;
  1137. X        xr = floor(active[j+1].x-.5);    /* right end of span */
  1138. X        if (xr>win->x1) xr = win->x1;
  1139. X        if (xl<=xr)
  1140. X        (*spanproc)(y, xl, xr);        /* draw pixels in span */
  1141. X        active[j].x += active[j].dx;    /* increment edge coords */
  1142. X        active[j+1].x += active[j+1].dx;
  1143. X    }
  1144. X    }
  1145. X}
  1146. X
  1147. Xstatic delete(i)        /* remove edge i from active list */
  1148. Xint i;
  1149. X{
  1150. X    int j;
  1151. X
  1152. X    for (j=0; j<nact && active[j].i!=i; j++);
  1153. X    if (j>=nact) return;    /* edge not in active list; happens at win->y0*/
  1154. X    nact--;
  1155. X    bcopy(&active[j+1], &active[j], (nact-j)*sizeof active[0]);
  1156. X}
  1157. X
  1158. Xstatic insert(i, y)        /* append edge i to end of active list */
  1159. Xint i, y;
  1160. X{
  1161. X    int j;
  1162. X    double dx;
  1163. X    Point2 *p, *q;
  1164. X
  1165. X    j = i<n-1 ? i+1 : 0;
  1166. X    if (pt[i].y < pt[j].y) {p = &pt[i]; q = &pt[j];}
  1167. X    else           {p = &pt[j]; q = &pt[i];}
  1168. X    /* initialize x position at intersection of edge with scanline y */
  1169. X    active[nact].dx = dx = (q->x-p->x)/(q->y-p->y);
  1170. X    active[nact].x = dx*(y+.5-p->y)+p->x;
  1171. X    active[nact].i = i;
  1172. X    nact++;
  1173. X}
  1174. X
  1175. X/* comparison routines for qsort */
  1176. Xcompare_ind(u, v) int *u, *v; {return pt[*u].y <= pt[*v].y ? -1 : 1;}
  1177. Xcompare_active(u, v) Edge *u, *v; {return u->x <= v->x ? -1 : 1;}
  1178. END_OF_FILE
  1179. if test 4931 -ne `wc -c <'ConcaveScan.c'`; then
  1180.     echo shar: \"'ConcaveScan.c'\" unpacked with wrong size!
  1181. fi
  1182. # end of 'ConcaveScan.c'
  1183. fi
  1184. if test -f 'Forms.c' -a "${1}" != "-c" ; then 
  1185.   echo shar: Will not clobber existing file \"'Forms.c'\"
  1186. else
  1187. echo shar: Extracting \"'Forms.c'\" \(5341 characters\)
  1188. sed "s/^X//" >'Forms.c' <<'END_OF_FILE'
  1189. X/*
  1190. XForms, Vectors,and Transforms
  1191. Xby Bob Wallis
  1192. Xfrom "Graphics Gems", Academic Press, 1990
  1193. X*/
  1194. X
  1195. X/*----------------------------------------------------------------------
  1196. XThe main program below is set up to solve the Bezier subdivision problem
  1197. Xin "Forms, Vectors, and Transforms".  The subroutines are useful in
  1198. Xsolving general problems which require manipulating matrices via exact
  1199. Xinteger arithmetic.  The intended application is validating or avoiding
  1200. Xtedious algebraic calculations.  As such, no thought was given to
  1201. Xefficiency.  
  1202. X----------------------------------------------------------------------*/
  1203. X#define ABS(x) ((x)>(0)? (x):(-x))
  1204. X#define N 4            /* size of matrices to deal with */
  1205. Xint     M[N][N] =        /* Bezier weights */
  1206. X{
  1207. X      1,   0,   0,   0,
  1208. X     -3,   3,   0,   0,
  1209. X      3,  -6,   3,   0,
  1210. X     -1,   3,  -3,   1,
  1211. X};
  1212. Xint     T[N][N] =    /* re-parameterization xform for top half */
  1213. X{
  1214. X      1,  -1,   1,  -1,
  1215. X      0,   2,  -4,   6,
  1216. X      0,   0,   4, -12,
  1217. X      0,   0,   0,   8
  1218. X};
  1219. Xmain ()
  1220. X{
  1221. X    int     i,
  1222. X            j,
  1223. X            scale,
  1224. X            gcd,
  1225. X            C[N][N],
  1226. X            S[N][N],
  1227. X            Madj[N][N],
  1228. X            Tadj[N][N],
  1229. X            Mdet,
  1230. X            Tdet;
  1231. X
  1232. X
  1233. X    Tdet = adjoint (T, Tadj);   /* inverse without division by */
  1234. X    Mdet = adjoint (M, Madj);   /* determinant of T and M */
  1235. X    matmult (Madj, Tadj, C);
  1236. X    matmult (C, M, S);        /* Madj*Tadj*M -> S */
  1237. X    scale = gcd = Mdet * Tdet;  /* scale factors of both determinants */
  1238. X    for (i = 0; i < N; i++)    /* find the greatest common */
  1239. X    {                /* demoninator of S and determinants */
  1240. X        for (j = 0; j < N; j++)
  1241. X            gcd = Gcd (gcd, S[i][j]);
  1242. X    }
  1243. X    scale /= gcd;        /* divide everything by gcd to get */
  1244. X    for (i = 0; i < N; i++)    /* matrix and scale factor in lowest */
  1245. X    {                /* integer terms possible */
  1246. X        for (j = 0; j < N; j++)
  1247. X            S[i][j] /= gcd;
  1248. X    }
  1249. X    printf ("scale factor = 1/%d  ", scale);
  1250. X    print_mat ("M=", M, N);     /* display the results */
  1251. X    print_mat ("T=", T, N);
  1252. X    print_mat ("S=", S, N);     /* subdivision matrix */
  1253. X    exit (0);
  1254. X}
  1255. XGcd (a, b)            /*returns greatest common demoninator */
  1256. Xint     a,            /* of (a,b) */
  1257. X           b;
  1258. X{
  1259. X    int        i,
  1260. X                r;
  1261. X    a = ABS (a);        /* force positive */
  1262. X    b = ABS (b);
  1263. X    if (a < b)            /* exchange so that a >= b */
  1264. X    {
  1265. X        i = b;
  1266. X        b = a;
  1267. X        a = i;
  1268. X    }
  1269. X    if     (b == 0)
  1270. X        return (a);    /* finished */
  1271. X    r = a % b;            /* remainder */
  1272. X    if (r == 0)
  1273. X        return (b);    /* finished */
  1274. X    else            /* recursive call */
  1275. X        return (Gcd (b, r));
  1276. X}
  1277. X
  1278. Xadjoint (A, Aadj)            /* returns determinant of A */
  1279. Xint     A[N][N],            /* input matrix */
  1280. X        Aadj[N][N];            /* output = adjoint of A */
  1281. X{                    /* must have N >= 3 */
  1282. X    int       i,
  1283. X            j,
  1284. X            I[N],            /* arrays of row and column indices */
  1285. X            J[N],
  1286. X            Isub[N],            /* sub-arrays of the above */
  1287. X            Jsub[N],
  1288. X            cofactor,
  1289. X            det;
  1290. X    if (N < 3)
  1291. X    {
  1292. X        printf ("must have N >= 3\n");
  1293. X        exit (1);
  1294. X    }
  1295. X    for (i = 0; i < N; i++)
  1296. X    {                    /* lookup tables to select a */
  1297. X        I[i] = i;        /* particular subset of */
  1298. X        J[i] = i;        /* rows and columns */
  1299. X    }
  1300. X    det = 0;
  1301. X    for (i = 0; i < N; i++)
  1302. X    {                    /* delete ith row */
  1303. X        subarray (I, Isub, N, i);
  1304. X        for (j = 0; j < N; j++)
  1305. X        {                /* delete jth column */
  1306. X                subarray (J, Jsub, N, j);
  1307. X                cofactor = determinant (A, Isub, Jsub, N - 1,(i + j) & 1);
  1308. X                if (j == 0)        /* use 0th column for det */
  1309. X                det += cofactor * A[i][0];
  1310. X            Aadj[j][i] = cofactor;
  1311. X        }
  1312. X       }
  1313. X    return (det);
  1314. X}
  1315. Xdeterminant (A, I, J, n, parity)/* actually gets a sub-determinant */
  1316. Xint     A[N][N],            /* input = entire matrix */
  1317. X           I[N],                /* row sub-array we want */
  1318. X        J[N],                /* col sub-array we want */
  1319. X        parity,                /* 1-> flip polarity */
  1320. X        n;                /* # elements in subarrays */
  1321. X
  1322. X{
  1323. X    int     i,
  1324. X            j,
  1325. X            det,
  1326. X            j_,
  1327. X            Jsub[N];
  1328. X    if (n <= 2)            /* call ourselves till we get down to */
  1329. X    {                /* a 2x2 matrix */
  1330. X        det =
  1331. X            (A[I[0]][J[0]] * A[I[1]][J[1]]) -
  1332. X            (A[I[1]][J[0]] * A[I[0]][J[1]]);
  1333. X        if (parity)
  1334. X        det = -det;
  1335. X        return (det);
  1336. X    }                    /* if (n <= 2) */
  1337. X    det = 0;                /* n > 2; call recursively */
  1338. X    i = I[0];                /* strike out 0th row */
  1339. X    for (j_ = 0; j_ < n; j_++)
  1340. X    {                    /* strike out jth column */
  1341. X        subarray (J, Jsub, n, j_);
  1342. X        j = J[j_];        /* I + 1 => struck out 0th row */
  1343. X        det += A[i][j] * determinant (A, I + 1, Jsub, n - 1, j_ & 1);
  1344. X    }
  1345. X    if (parity)
  1346. X        det = -det;
  1347. X    return (det);
  1348. X}
  1349. Xsubarray (src, dest, n, k)        /* strike out kth row/column */
  1350. Xint      *src,                /*     source array of n indices */
  1351. X         *dest,                /* dest array formed by deleting k  */
  1352. X            n,
  1353. X            k;
  1354. X{
  1355. X    int     i;
  1356. X    for (i = 0; i < n; i++, src++)
  1357. X        if (i != k)            /* skip over k */
  1358. X            *dest++ = *src;
  1359. X}
  1360. X
  1361. Xmatmult (A, B, C)                /* C = A*B */
  1362. Xint     A[N][N],
  1363. X        B[N][N],
  1364. X        C[N][N];
  1365. X{
  1366. X    int        i,
  1367. X        j,
  1368. X        k,
  1369. X        sum;
  1370. X    for (i = 0; i < N; i++)
  1371. X    {
  1372. X        for (k = 0; k < N; k++)
  1373. X        {
  1374. X            sum = 0;
  1375. X            for (j = 0; j < N; j++)
  1376. X                sum += A[i][j] * B[j][k];
  1377. X            C[i][k] = sum;
  1378. X        }
  1379. X    }
  1380. X}
  1381. Xprint_mat (string, mat, n)
  1382. Xchar   *string;
  1383. Xint     mat[N][N],
  1384. X         n;
  1385. X{
  1386. X    int     i,
  1387. X             j;
  1388. X    printf ("%s\n", string);
  1389. X    for (i = 0; i < n; i++)
  1390. X    {
  1391. X        for (j = 0; j < n; j++)
  1392. X        printf (" %8ld", mat[i][j]);
  1393. X    printf ("\n");
  1394. X    }
  1395. X}
  1396. END_OF_FILE
  1397. if test 5341 -ne `wc -c <'Forms.c'`; then
  1398.     echo shar: \"'Forms.c'\" unpacked with wrong size!
  1399. fi
  1400. # end of 'Forms.c'
  1401. fi
  1402. if test -f 'Interleave.c' -a "${1}" != "-c" ; then 
  1403.   echo shar: Will not clobber existing file \"'Interleave.c'\"
  1404. else
  1405. echo shar: Extracting \"'Interleave.c'\" \(6497 characters\)
  1406. sed "s/^X//" >'Interleave.c' <<'END_OF_FILE'
  1407. X/*
  1408. XBit Interleaving for Quad- or Octrees
  1409. Xby Clifford A. Shaffer
  1410. Xfrom "Graphics Gems", Academic Press, 1990
  1411. X*/
  1412. X
  1413. X#include "GraphicsGems.h"
  1414. X#define B_MAX_DEPTH 14    /* maximum depth allowed */
  1415. X
  1416. X/* byteval is the lookup table for coordinate interleaving.  Given a
  1417. X   4 bit portion of the (x, y) coordinates, return the bit interleaving.
  1418. X   Notice that this table looks like the order in which the pixels of
  1419. X   a 16 X 16 pixel image would be visited. */
  1420. Xint byteval[16][16] = 
  1421. X    {  0,  1,  4,  5, 16, 17, 20, 21, 64, 65, 68, 69, 80, 81, 84, 85,
  1422. X       2,  3,  6,  7, 18, 19, 22, 23, 66, 67, 70, 71, 82, 83, 86, 87,
  1423. X       8,  9, 12, 13, 24, 25, 28, 29, 72, 73, 76, 77, 88, 89, 92, 93,
  1424. X      10, 11, 14, 15, 26, 27, 30, 31, 74, 75, 78, 79, 90, 91, 94, 95,
  1425. X      32, 33, 36, 37, 48, 49, 52, 53, 96, 97,100,101,112,113,116,117,
  1426. X      34, 35, 38, 39, 50, 51, 54, 55, 98, 99,102,103,114,115,118,119,
  1427. X      40, 41, 44, 45, 56, 57, 60, 61,104,105,108,109,120,121,124,125,
  1428. X      42, 43, 46, 47, 58, 59, 62, 63,106,107,110,111,122,123,126,127,
  1429. X     128,129,132,133,144,145,148,149,192,193,196,197,208,209,212,213,
  1430. X     130,131,134,135,146,147,150,151,194,195,198,199,210,211,214,215,
  1431. X     136,137,140,141,152,153,156,157,200,201,204,205,216,217,220,221,
  1432. X     138,139,142,143,154,155,158,159,202,203,206,207,218,219,222,223,
  1433. X     160,161,164,165,176,177,180,181,224,225,228,229,240,241,244,245,
  1434. X     162,163,166,167,178,179,182,183,226,227,230,231,242,243,246,247,
  1435. X     168,169,172,173,184,185,188,189,232,233,236,237,248,249,252,253,
  1436. X     170,171,174,175,186,187,190,191,234,235,238,239,250,251,254,255};
  1437. X
  1438. X/* bytemask is the mask for byte interleaving - masks out the 
  1439. X   non-significant bit positions.  This is determined by the
  1440. X   depth of the node. For example, a node of depth 0 is at the root.
  1441. X  Thus, there are no branchs and no bits are significant. 
  1442. X  The bottom 4 bits (the depth) are always retained. 
  1443. X  Values are in octal notation. */
  1444. Xint bytemask[B_MAX_DEPTH + 1] = {017,
  1445. X     030000000017,036000000017,037400000017,037700000017,
  1446. X     037760000017,037774000017,037777000017,037777600017,
  1447. X     037777740017,037777770017,037777776017,037777777417,
  1448. X     037777777717,037777777777};
  1449. X
  1450. X
  1451. Xlong *interleave(addr, x, y, depth, max_depth)
  1452. X/* Return the interleaved code for a quadtree node at depth depth 
  1453. Xwhose upper left hand corner has coordinates (x, y) in a tree with maximum
  1454. Xdepth max_depth.  This function receives and returns a 
  1455. Xpointer to addr, which is either a long interger or (more typically)
  1456. Xan array of long integers whose first integer contains the 
  1457. Xinterleaved code. */
  1458. X long *addr;
  1459. X int max_depth, depth;
  1460. X int x, y;
  1461. X{
  1462. X
  1463. X/* Scale x, y values to be consistent with maximum coord size */
  1464. X/* and depth of tree */
  1465. X x <<= (B_MAX_DEPTH - max_depth);
  1466. X y <<= (B_MAX_DEPTH - max_depth);
  1467. X
  1468. X/* calculate the bit interleaving of the x, y values that have now
  1469. X   been appropriately shifted, and place this interleave in the address
  1470. X   portion of addr.  Note that the binary representations of x and y are
  1471. X   being processed from right to left */
  1472. X
  1473. X *addr = depth;
  1474. X *addr |= byteval[y & 03][x & 03] << 4;
  1475. X *addr |= byteval[(y >> 2) & 017][(x >> 2) & 017] << 8;
  1476. X *addr |= byteval[(y >> 6) & 017][(x >> 6) & 017] << 16;
  1477. X *addr |= byteval[(y >> 10) & 017][(x >> 10) & 017] << 24;
  1478. X *addr &= bytemask[depth];
  1479. X
  1480. X/* if there were unused portions of the x and y addresses then  */
  1481. X/* the address was too large for the depth values given.  */
  1482. X/*  Return address built */
  1483. X return (addr);
  1484. X}
  1485. X
  1486. X
  1487. X
  1488. X/* The next two arrays are used in calculating the (x, y) coordinates
  1489. X   of the upper left-hand corner of a node from its bit interleaved
  1490. X   address.  Given an 8 bit number, the arrays return the effect of
  1491. X   removing every other bit (the y bits preceed the x bits). */
  1492. X
  1493. Xint xval[256] = { 0, 1, 0, 1, 2, 3, 2, 3, 0, 1, 0, 1, 2, 3, 2, 3,
  1494. X                 4, 5, 4, 5, 6, 7, 6, 7, 4, 5, 4, 5, 6, 7, 6, 7,
  1495. X                 0, 1, 0, 1, 2, 3, 2, 3, 0, 1, 0, 1, 2, 3, 2, 3,
  1496. X                 4, 5, 4, 5, 6, 7, 6, 7, 4, 5, 4, 5, 6, 7, 6, 7,
  1497. X                 8, 9, 8, 9,10,11,10,11, 8, 9, 8, 9,10,11,10,11,
  1498. X                12,13,12,13,14,15,14,15,12,13,12,13,14,15,14,15,
  1499. X                 8, 9, 8, 9,10,11,10,11, 8, 9, 8, 9,10,11,10,11,
  1500. X                12,13,12,13,14,15,14,15,12,13,12,13,14,15,14,15,
  1501. X                 0, 1, 0, 1, 2, 3, 2, 3, 0, 1, 0, 1, 2, 3, 2, 3,
  1502. X                 4, 5, 4, 5, 6, 7, 6, 7, 4, 5, 4, 5, 6, 7, 6, 7,
  1503. X                 0, 1, 0, 1, 2, 3, 2, 3, 0, 1, 0, 1, 2, 3, 2, 3,
  1504. X                 4, 5, 4, 5, 6, 7, 6, 7, 4, 5, 4, 5, 6, 7, 6, 7,
  1505. X                 8, 9, 8, 9,10,11,10,11, 8, 9, 8, 9,10,11,10,11,
  1506. X                12,13,12,13,14,15,14,15,12,13,12,13,14,15,14,15,
  1507. X                 8, 9, 8, 9,10,11,10,11, 8, 9, 8, 9,10,11,10,11,
  1508. X                12,13,12,13,14,15,14,15,12,13,12,13,14,15,14,15};
  1509. X
  1510. X
  1511. Xint yval[256] = { 0, 0, 1, 1, 0, 0, 1, 1, 2, 2, 3, 3, 2, 2, 3, 3,
  1512. X              0, 0, 1, 1, 0, 0, 1, 1, 2, 2, 3, 3, 2, 2, 3, 3,
  1513. X              4, 4, 5, 5, 4, 4, 5, 5, 6, 6, 7, 7, 6, 6, 7, 7,
  1514. X              4, 4, 5, 5, 4, 4, 5, 5, 6, 6, 7, 7, 6, 6, 7, 7,
  1515. X              0, 0, 1, 1, 0, 0, 1, 1, 2, 2, 3, 3, 2, 2, 3, 3,
  1516. X              0, 0, 1, 1, 0, 0, 1, 1, 2, 2, 3, 3, 2, 2, 3, 3,
  1517. X              4, 4, 5, 5, 4, 4, 5, 5, 6, 6, 7, 7, 6, 6, 7, 7,
  1518. X              4, 4, 5, 5, 4, 4, 5, 5, 6, 6, 7, 7, 6, 6, 7, 7,
  1519. X              8, 8, 9, 9, 8, 8, 9, 9,10,10,11,11,10,10,11,11,
  1520. X              8, 8, 9, 9, 8, 8, 9, 9,10,10,11,11,10,10,11,11,
  1521. X             12,12,13,13,12,12,13,13,14,14,15,15,14,14,15,15,
  1522. X             12,12,13,13,12,12,13,13,14,14,15,15,14,14,15,15,
  1523. X              8, 8, 9, 9, 8, 8, 9, 9,10,10,11,11,10,10,11,11,
  1524. X              8, 8, 9, 9, 8, 8, 9, 9,10,10,11,11,10,10,11,11,
  1525. X             12,12,13,13,12,12,13,13,14,14,15,15,14,14,15,15,
  1526. X             12,12,13,13,12,12,13,13,14,14,15,15,14,14,15,15};
  1527. X
  1528. X
  1529. X
  1530. X
  1531. X
  1532. Xint getx(addr, max_depth)
  1533. X/* Return the x coordinate of the upper left hand corner of addr for a
  1534. X   tree with maximum depth max_depth. */
  1535. X long *addr;
  1536. X int max_depth;
  1537. X{
  1538. X register x;
  1539. X
  1540. X x = xval[(*addr >> 4) & 017];
  1541. X x |= xval[(*addr >> 8) & 0377] << 2;
  1542. X x |= xval[(*addr >> 16) & 0377] << 6;
  1543. X x |= xval[(*addr >> 24) & 0377] << 10;
  1544. X x >>= B_MAX_DEPTH - max_depth;
  1545. X return (x);
  1546. X}
  1547. X
  1548. X
  1549. X
  1550. Xint QKy(addr, max_depth)
  1551. X/* Return the y coordinate of the upper left hand corner of addr for a
  1552. X   tree with maximum depth max_depth. */
  1553. X
  1554. X long *addr;
  1555. X int max_depth;
  1556. X{
  1557. X register y;
  1558. X
  1559. X y = yval[(*addr >> 4) & 017];
  1560. X y |= yval[(*addr >> 8) & 0377] << 2;
  1561. X y |= yval[(*addr >> 16) & 0377] << 6;
  1562. X y |= yval[(*addr >> 24) & 0377] << 10;
  1563. X y >>= B_MAX_DEPTH - max_depth;
  1564. X return (y);
  1565. X}
  1566. X
  1567. Xint getdepth(addr)
  1568. X/* Return the depth of the node.  Simply return the bottom 4 bits. */
  1569. X
  1570. X long *addr;
  1571. X{
  1572. X
  1573. X return(*addr & 017);
  1574. X}
  1575. END_OF_FILE
  1576. if test 6497 -ne `wc -c <'Interleave.c'`; then
  1577.     echo shar: \"'Interleave.c'\" unpacked with wrong size!
  1578. fi
  1579. # end of 'Interleave.c'
  1580. fi
  1581. if test -f 'Sturm/sturm.c' -a "${1}" != "-c" ; then 
  1582.   echo shar: Will not clobber existing file \"'Sturm/sturm.c'\"
  1583. else
  1584. echo shar: Extracting \"'Sturm/sturm.c'\" \(5198 characters\)
  1585. sed "s/^X//" >'Sturm/sturm.c' <<'END_OF_FILE'
  1586. X
  1587. X/*
  1588. X * sturm.c
  1589. X *
  1590. X *    the functions to build and evaluate the Sturm sequence
  1591. X */
  1592. X#include <math.h>
  1593. X#include <stdio.h>
  1594. X#include "solve.h"
  1595. X
  1596. X/*
  1597. X * modp
  1598. X *
  1599. X *    calculates the modulus of u(x) / v(x) leaving it in r, it
  1600. X *  returns 0 if r(x) is a constant.
  1601. X *  note: this function assumes the leading coefficient of v 
  1602. X *    is 1 or -1
  1603. X */
  1604. Xstatic int
  1605. Xmodp(u, v, r)
  1606. X    poly    *u, *v, *r;
  1607. X{
  1608. X    int        k, j;
  1609. X    double    *nr, *end, *uc;
  1610. X
  1611. X    nr = r->coef;
  1612. X    end = &u->coef[u->ord];
  1613. X
  1614. X    uc = u->coef;
  1615. X    while (uc <= end)
  1616. X            *nr++ = *uc++;
  1617. X
  1618. X    if (v->coef[v->ord] < 0.0) {
  1619. X
  1620. X
  1621. X            for (k = u->ord - v->ord - 1; k >= 0; k -= 2)
  1622. X                r->coef[k] = -r->coef[k];
  1623. X
  1624. X            for (k = u->ord - v->ord; k >= 0; k--)
  1625. X                for (j = v->ord + k - 1; j >= k; j--)
  1626. X                    r->coef[j] = -r->coef[j] - r->coef[v->ord + k]
  1627. X                    * v->coef[j - k];
  1628. X    } else {
  1629. X            for (k = u->ord - v->ord; k >= 0; k--)
  1630. X                for (j = v->ord + k - 1; j >= k; j--)
  1631. X                r->coef[j] -= r->coef[v->ord + k] * v->coef[j - k];
  1632. X    }
  1633. X
  1634. X    k = v->ord - 1;
  1635. X    while (k >= 0 && fabs(r->coef[k]) < SMALL_ENOUGH) {
  1636. X        r->coef[k] = 0.0;
  1637. X        k--;
  1638. X    }
  1639. X
  1640. X    r->ord = (k < 0) ? 0 : k;
  1641. X
  1642. X    return(r->ord);
  1643. X}
  1644. X
  1645. X/*
  1646. X * buildsturm
  1647. X *
  1648. X *    build up a sturm sequence for a polynomial in smat, returning
  1649. X * the number of polynomials in the sequence
  1650. X */
  1651. Xint
  1652. Xbuildsturm(ord, sseq)
  1653. X    int        ord;
  1654. X    poly    *sseq;
  1655. X{
  1656. X    int        i;
  1657. X    double    f, *fp, *fc;
  1658. X    poly    *sp;
  1659. X
  1660. X    sseq[0].ord = ord;
  1661. X    sseq[1].ord = ord - 1;
  1662. X
  1663. X
  1664. X    /*
  1665. X     * calculate the derivative and normalise the leading
  1666. X     * coefficient.
  1667. X     */
  1668. X    f = fabs(sseq[0].coef[ord] * ord);
  1669. X    fp = sseq[1].coef;
  1670. X    fc = sseq[0].coef + 1;
  1671. X    for (i = 1; i <= ord; i++)
  1672. X            *fp++ = *fc++ * i / f;
  1673. X
  1674. X    /*
  1675. X     * construct the rest of the Sturm sequence
  1676. X     */
  1677. X    for (sp = sseq + 2; modp(sp - 2, sp - 1, sp); sp++) {
  1678. X
  1679. X        /*
  1680. X         * reverse the sign and normalise
  1681. X         */
  1682. X        f = -fabs(sp->coef[sp->ord]);
  1683. X        for (fp = &sp->coef[sp->ord]; fp >= sp->coef; fp--)
  1684. X                *fp /= f;
  1685. X    }
  1686. X
  1687. X    sp->coef[0] = -sp->coef[0];    /* reverse the sign */
  1688. X
  1689. X    return(sp - sseq);
  1690. X}
  1691. X
  1692. X/*
  1693. X * numroots
  1694. X *
  1695. X *    return the number of distinct real roots of the polynomial
  1696. X * described in sseq.
  1697. X */
  1698. Xint
  1699. Xnumroots(np, sseq, atneg, atpos)
  1700. X        int        np;
  1701. X        poly    *sseq;
  1702. X        int        *atneg, *atpos;
  1703. X{
  1704. X        int        atposinf, atneginf;
  1705. X        poly    *s;
  1706. X        double    f, lf;
  1707. X
  1708. X        atposinf = atneginf = 0;
  1709. X
  1710. X
  1711. X    /*
  1712. X     * changes at positve infinity
  1713. X     */
  1714. X    lf = sseq[0].coef[sseq[0].ord];
  1715. X
  1716. X    for (s = sseq + 1; s <= sseq + np; s++) {
  1717. X            f = s->coef[s->ord];
  1718. X            if (lf == 0.0 || lf * f < 0)
  1719. X                atposinf++;
  1720. X        lf = f;
  1721. X    }
  1722. X
  1723. X    /*
  1724. X     * changes at negative infinity
  1725. X     */
  1726. X    if (sseq[0].ord & 1)
  1727. X            lf = -sseq[0].coef[sseq[0].ord];
  1728. X    else
  1729. X            lf = sseq[0].coef[sseq[0].ord];
  1730. X
  1731. X    for (s = sseq + 1; s <= sseq + np; s++) {
  1732. X            if (s->ord & 1)
  1733. X                f = -s->coef[s->ord];
  1734. X            else
  1735. X                f = s->coef[s->ord];
  1736. X            if (lf == 0.0 || lf * f < 0)
  1737. X                atneginf++;
  1738. X            lf = f;
  1739. X    }
  1740. X
  1741. X    *atneg = atneginf;
  1742. X    *atpos = atposinf;
  1743. X
  1744. X    return(atneginf - atposinf);
  1745. X}
  1746. X
  1747. X/*
  1748. X * numchanges
  1749. X *
  1750. X *    return the number of sign changes in the Sturm sequence in
  1751. X * sseq at the value a.
  1752. X */
  1753. Xint
  1754. Xnumchanges(np, sseq, a)
  1755. X    int        np;
  1756. X    poly    *sseq;
  1757. X    double    a;
  1758. X
  1759. X{
  1760. X    int        changes;
  1761. X    double    f, lf;
  1762. X    poly    *s;
  1763. X
  1764. X    changes = 0;
  1765. X
  1766. X    lf = evalpoly(sseq[0].ord, sseq[0].coef, a);
  1767. X
  1768. X    for (s = sseq + 1; s <= sseq + np; s++) {
  1769. X            f = evalpoly(s->ord, s->coef, a);
  1770. X            if (lf == 0.0 || lf * f < 0)
  1771. X                changes++;
  1772. X            lf = f;
  1773. X    }
  1774. X
  1775. X    return(changes);
  1776. X}
  1777. X
  1778. X/*
  1779. X * sbisect
  1780. X *
  1781. X *    uses a bisection based on the sturm sequence for the polynomial
  1782. X * described in sseq to isolate intervals in which roots occur,
  1783. X * the roots are returned in the roots array in order of magnitude.
  1784. X */
  1785. Xsbisect(np, sseq, min, max, atmin, atmax, roots)
  1786. X    int        np;
  1787. X    poly    *sseq;
  1788. X    double    min, max;
  1789. X    int        atmin, atmax;
  1790. X    double    *roots;
  1791. X{
  1792. X    double    mid;
  1793. X    int        n1, n2, its, atmid, nroot;
  1794. X
  1795. X    if ((nroot = atmin - atmax) == 1) {
  1796. X
  1797. X        /*
  1798. X         * first try a less expensive technique.
  1799. X         */
  1800. X        if (modrf(sseq->ord, sseq->coef, min, max, &roots[0]))
  1801. X            return;
  1802. X
  1803. X
  1804. X        /*
  1805. X         * if we get here we have to evaluate the root the hard
  1806. X         * way by using the Sturm sequence.
  1807. X         */
  1808. X        for (its = 0; its < MAXIT; its++) {
  1809. X                mid = (min + max) / 2;
  1810. X
  1811. X                atmid = numchanges(np, sseq, mid);
  1812. X
  1813. X                if (fabs(mid) > RELERROR) {
  1814. X                    if (fabs((max - min) / mid) < RELERROR) {
  1815. X                        roots[0] = mid;
  1816. X                        return;
  1817. X                    }
  1818. X                } else if (fabs(max - min) < RELERROR) {
  1819. X                    roots[0] = mid;
  1820. X                    return;
  1821. X                }
  1822. X
  1823. X                if ((atmin - atmid) == 0)
  1824. X                    min = mid;
  1825. X                else
  1826. X                    max = mid;
  1827. X            }
  1828. X
  1829. X        if (its == MAXIT) {
  1830. X                fprintf(stderr, "sbisect: overflow min %f max %f\
  1831. X                    diff %e nroot %d n1 %d n2 %d\n",
  1832. X                    min, max, max - min, nroot, n1, n2);
  1833. X            roots[0] = mid;
  1834. X        }
  1835. X
  1836. X        return;
  1837. X    }
  1838. X
  1839. X    /*
  1840. X     * more than one root in the interval, we have to bisect...
  1841. X     */
  1842. X    for (its = 0; its < MAXIT; its++) {
  1843. X
  1844. X            mid = (min + max) / 2;
  1845. X
  1846. X            atmid = numchanges(np, sseq, mid);
  1847. X
  1848. X            n1 = atmin - atmid;
  1849. X            n2 = atmid - atmax;
  1850. X
  1851. X
  1852. X            if (n1 != 0 && n2 != 0) {
  1853. X                sbisect(np, sseq, min, mid, atmin, atmid, roots);
  1854. X                sbisect(np, sseq, mid, max, atmid, atmax, &roots[n1]);
  1855. X                break;
  1856. X            }
  1857. X
  1858. X            if (n1 == 0)
  1859. X                min = mid;
  1860. X            else
  1861. X                max = mid;
  1862. X    }
  1863. X
  1864. X    if (its == MAXIT) {
  1865. X            fprintf(stderr, "sbisect: roots too close together\n");
  1866. X            fprintf(stderr, "sbisect: overflow min %f max %f diff %e\
  1867. X                nroot %d n1 %d n2 %d\n",
  1868. X                min, max, max - min, nroot, n1, n2);
  1869. X            for (n1 = atmax; n1 < atmin; n1++)
  1870. X            roots[n1 - atmax] = mid;
  1871. X    }
  1872. X}
  1873. END_OF_FILE
  1874. if test 5198 -ne `wc -c <'Sturm/sturm.c'`; then
  1875.     echo shar: \"'Sturm/sturm.c'\" unpacked with wrong size!
  1876. fi
  1877. # end of 'Sturm/sturm.c'
  1878. fi
  1879. echo shar: End of archive 4 \(of 5\).
  1880. cp /dev/null ark4isdone
  1881. MISSING=""
  1882. for I in 1 2 3 4 5 ; do
  1883.     if test ! -f ark${I}isdone ; then
  1884.     MISSING="${MISSING} ${I}"
  1885.     fi
  1886. done
  1887. if test "${MISSING}" = "" ; then
  1888.     echo You have unpacked all 5 archives.
  1889.     rm -f ark[1-9]isdone
  1890. else
  1891.     echo You still need to unpack the following archives:
  1892.     echo "        " ${MISSING}
  1893. fi
  1894. ##  End of shell archive.
  1895. exit 0
  1896.  
  1897.